| python求解线性规划问题 | 您所在的位置：网站首页 › 豆汁 纳豆 › python求解线性规划问题 |
python求解线性规划问题
|
两阶段法解一般的标准形式的线性规划问题
1、两阶段法是什么
单纯形法的关键在与如何找到初始基可行解，而两阶段法通过添加人工变量构造第一阶段的另一个线性规划问题，使得基可行解易求得，然后通过迭代得到第一阶段的最优解，再讨论原问题的最优解情况。 2、例题原LP：（已经化为标准形式） m a x z = x 1 + 3 x 2 − x 3 s . t . { x 1 + x 2 + 2 x 3 + x 4 = 4 − x 1 + 2 x 2 + x 3 + x 4 = 4 3 x 1 + 3 x 3 + x 4 = 4 x i ≥ 0 i = 1 , 2 , 3 , 4 max \ \ z=x_1+3x_2-x_3 \\ s.t.\begin{cases} x_1+x_2+2x_3+x_4=4 \\ -x_1+2x_2+x_3+x_4=4 \\ 3x_1+3x_3+x_4=4 \\ x_i \geq0 \ \ \ \ \ \ \ i=1,2,3,4\end{cases} max z=x1+3x2−x3s.t.⎩⎪⎪⎪⎨⎪⎪⎪⎧x1+x2+2x3+x4=4−x1+2x2+x3+x4=43x1+3x3+x4=4xi≥0 i=1,2,3,4 LP1: m a x z 1 = − ( x 5 + x 6 + x 7 ) s . t . { x 1 + x 2 + 2 x 3 + x 4 + x 5 = 4 − x 1 + 2 x 2 + x 3 + x 4 + x 6 = 4 3 x 1 + 3 x 3 + x 4 + x 7 = 4 x i ≥ 0 i = 1 , 2 , 3 , 4 , 5 , 6 , 7 max \ \ z_1=-(x_5+x_6+x_7) \\ s.t.\begin{cases} x_1+x_2+2x_3+x_4+x_5=4 \\ -x_1+2x_2+x_3+x_4+x_6=4 \\ 3x_1+3x_3+x_4+x_7=4 \\ x_i \geq0 \ \ \ \ \ \ \ i=1,2,3,4,5,6,7\end{cases} max z1=−(x5+x6+x7)s.t.⎩⎪⎪⎪⎨⎪⎪⎪⎧x1+x2+2x3+x4+x5=4−x1+2x2+x3+x4+x6=43x1+3x3+x4+x7=4xi≥0 i=1,2,3,4,5,6,7 易知LP1若有最优解，且最优解中 X M = ( x 5 , x 6 , x 7 ) = ( 0 , 0 , 0 ) X_M=(x_5,x_6,x_7)=(0,0,0) XM=(x5,x6,x7)=(0,0,0)那么原问题就有基本可行解。 3、python实现 a、获取输入 def getinput(): global m,n m = int(input('输入约束的个数'))# m个约束条件，n个变量 n = int(input('输入变量的个数')) a = [] #约束系数矩阵 for i in range(m): a.append( [eval(input('a({}{})='.format(i+1, j+1))) for j in range(n)] ) e = [list(i) for i in np.diag([1]*m)] #人工变量系数矩阵（m阶单位阵) b = [eval(input('b({})='.format(i+1))) for i in range(m)] #约束条件右端常数 r_1 = [] #LP_1检验数 for i in range(n): r_1.append(sum([a[j][i] for j in range(m)])) r_1 += [0]*m r_1.append(sum(b)) r = [eval(input('r({})='.format(i+1))) for i in range(n)]+[0]*(m+1) vect = [n+i+1 for i in range(m)] # 基变量足标 a = [i+j+[k] for i, j, k in zip(a, e, b)]+[r_1]+[r] return a,vect b、判断是否需要转轴 def judge(matrix): if max(matrix[-2][:-1]) } #用来寻找满足min的变量 for i in a[:-2]: if i[index] >0: d[i[-1]/i[index]] = a.index(i) pivot = d[min(d)] matrix[pivot] = matrix[pivot]/matrix[pivot][index] for i in range(m+2): if i != pivot: matrix[i] = matrix[i] - matrix[i][index]*matrix[pivot] vect[pivot] = index+1 #基变量的变动 a = [list(i) for i in matrix] def pr(matrix,vect): #输出单纯形表 print('*'*20) #单纯为了好看 print('XB',end='\t') for i in range(n+m): print('X_{}'.format(i+1),end='\t') print('b') for i in range(m+2): if i m: print('(x_{})'.format(vect[i]),end='\t') else: print('x_{}'.format(vect[i]),end='\t') elif i == m: print('r_1',end='\t') elif i == m+1: print('r',end='\t') for j in matrix[i]: print(f(str(j)).limit_denominator(),end='\t') #输出分数形式 print(end='\n') def print_solution(matrix): print('*'*20) for i in range(n+m): print('x{}*={}'.format(i+1,matrix[-2][i]),end='，') print('\nz*={}'.format(-matrix[-2][-1])) def main(): global a,matrix,vect a,vect = getinput() matrix = np.array(a, dtype=np.float64) #array可以方便地进行整行的操作，而列表可以方便索引（其实就是我不会array的索引23333） pr(matrix,vect) while judge(matrix): trans() pr(matrix,vect) print_solution(matrix) main()输入输出示例： |
【本文地址】